Near Analysis Null Models



What Was Tested: Do crop yields increase from randomly placed GCOs globally? (T# indicated random GCO)



What This Allows Me To Do: Identify if yield-distance relationships (near analysis) are a function of evolutionary escape or other variables.



What Are The Model Assumptions: GLS models are used as they have less assumptions about error structure and heteroscadicity.

#Barley -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
barley.null<-null %>% select(mean_barle,rescale_ND)
barley.null<-barley.null %>% filter(mean_barle > 0, na.rm=TRUE)
barley.null$logBarley<-log10(barley.null$mean_barle)

#Plot
ggplot(barley.null, aes(x=rescale_ND, y=logBarley)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Barley Near Null Model, T01

Barley Near Null Model, T01 

###  Model - Barley T01
# regular OLS no variance structure
barley.ols <- gls(logBarley ~ rescale_ND, data = barley.null)
# varFixed (variance changes linearly with X)
barley.fixed <- update(barley.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
barley.power <- update(barley.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
barley.exp <- update(barley.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
barley.ConstPower <- update(barley.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(barley.ols, barley.fixed, barley.power, barley.ConstPower, barley.exp)
##                   df      AIC
## barley.ols         3 19860.72
## barley.fixed       3 20649.60
## barley.power       4 19818.06
## barley.ConstPower  5 19818.25
## barley.exp         4 19818.39
#varPower
summary(barley.power)
## Generalized least squares fit by REML
##   Model: logBarley ~ rescale_ND 
##   Data: barley.null 
##        AIC      BIC    logLik
##   19818.06 19847.35 -9905.029
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##    power 
## 0.101525 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept)  0.09030799 0.016095330  5.610820       0
## rescale_ND  -0.01134541 0.001242831 -9.128682       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.939
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.2685673 -0.4483371  0.2915946  0.7273243  1.5829412 
## 
## Residual standard error: 0.4567186 
## Degrees of freedom: 11197 total; 11195 residual

#Cassava -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
cassava.null<-null %>% select(mean_cassa,rescale_ND)
cassava.null<-cassava.null %>% filter(mean_cassa > 0, na.rm=TRUE)
cassava.null$logCassava<-log10(cassava.null$mean_cassa)

#Plot
ggplot(cassava.null, aes(x=rescale_ND, y=logCassava)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Cassava Near Null Model, T01

Cassava Near Null Model, T01 

###  Model - Cassava T01
# regular OLS no variance structure
cassava.ols <- gls(logCassava ~ rescale_ND, data = cassava.null)
# varFixed (variance changes linearly with X)
cassava.fixed <- update(cassava.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
cassava.power <- update(cassava.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
cassava.exp <- update(cassava.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
cassava.ConstPower <- update(cassava.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(cassava.ols, cassava.fixed, cassava.power, cassava.ConstPower, cassava.exp)
##                    df       AIC
## cassava.ols         3  9128.050
## cassava.fixed       3 10246.056
## cassava.power       4  9129.169
## cassava.ConstPower  5  9131.167
## cassava.exp         4  9127.686
#varExp
summary(cassava.exp)
## Generalized least squares fit by REML
##   Model: logCassava ~ rescale_ND 
##   Data: cassava.null 
##        AIC      BIC    logLik
##   9127.686 9154.688 -4559.843
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       expon 
## 0.002698391 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept) 0.6418420 0.015944886 40.25379       0
## rescale_ND  0.0080824 0.001226775  6.58837       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.92 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -5.3384659 -0.3436627  0.2918642  0.7277743  1.5659812 
## 
## Residual standard error: 0.4814093 
## Degrees of freedom: 6316 total; 6314 residual

#Groundnut -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
groundnut.null<-null %>% select(mean_groun,rescale_ND)
groundnut.null<-groundnut.null %>% filter(mean_groun > 0, na.rm=TRUE)
groundnut.null$logGroundnut<-log10(groundnut.null$mean_groun)

#Plot
ggplot(groundnut.null, aes(x=rescale_ND, y=logGroundnut)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Groundnut Near Null Model, T01

Groundnut Near Null Model, T01 

###  Model - Groundnut T01
# regular OLS no variance structure
groundnut.ols <- gls(logGroundnut ~ rescale_ND, data = groundnut.null)
# varFixed (variance changes linearly with X)
groundnut.fixed <- update(groundnut.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
groundnut.power <- update(groundnut.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
groundnut.exp <- update(groundnut.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
groundnut.ConstPower <- update(groundnut.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(groundnut.ols, groundnut.fixed, groundnut.power, groundnut.ConstPower, groundnut.exp)
##                      df      AIC
## groundnut.ols         3 13518.83
## groundnut.fixed       3 14568.52
## groundnut.power       4 13516.17
## groundnut.ConstPower  5 13518.17
## groundnut.exp         4 13520.29
#varPower
summary(groundnut.power)
## Generalized least squares fit by REML
##   Model: logGroundnut ~ rescale_ND 
##   Data: groundnut.null 
##        AIC      BIC    logLik
##   13516.17 13544.27 -6754.084
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##     power 
## 0.0340043 
## 
## Coefficients:
##                  Value   Std.Error    t-value p-value
## (Intercept) -0.2628995 0.016587741 -15.849027       0
## rescale_ND   0.0050640 0.001237187   4.093157       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.933
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.4907423 -0.4119282  0.2553681  0.6928837  1.6701106 
## 
## Residual standard error: 0.5010125 
## Degrees of freedom: 8322 total; 8320 residual

#Maize -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
maize.null<-null %>% select(mean_maize,rescale_ND)
maize.null<-maize.null %>% filter(mean_maize > 0, na.rm=TRUE)
maize.null$logMaize<-log10(maize.null$mean_maize)

#Plot
ggplot(maize.null, aes(x=rescale_ND, y=logMaize)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Maize Near Null Model, T01

Maize Near Null Model, T01 

###  Model - Maize T01
# regular OLS no variance structure
maize.ols <- gls(logMaize ~ rescale_ND, data = maize.null)
# varFixed (variance changes linearly with X)
maize.fixed <- update(maize.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
maize.power <- update(maize.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x) - Does not Converge
#maize.exp <- update(maize.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
maize.ConstPower <- update(maize.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(maize.ols, maize.fixed, maize.power, maize.ConstPower) #,maize.exp)
##                  df      AIC
## maize.ols         3 24154.27
## maize.fixed       3 26364.94
## maize.power       4 24151.90
## maize.ConstPower  5 24153.91
#varPower
summary(maize.power)
## Generalized least squares fit by REML
##   Model: logMaize ~ rescale_ND 
##   Data: maize.null 
##       AIC      BIC    logLik
##   24151.9 24182.03 -12071.95
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       power 
## -0.02496171 
## 
## Coefficients:
##                  Value   Std.Error   t-value p-value
## (Intercept)  0.3224096 0.013719180  23.50064       0
## rescale_ND  -0.0156234 0.001062368 -14.70622       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.933
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.6035233 -0.4830642  0.2177102  0.6990680  2.0207779 
## 
## Residual standard error: 0.6155705 
## Degrees of freedom: 13793 total; 13791 residual

#Millet -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
millet.null<-null %>% select(mean_mille,rescale_ND)
millet.null<-millet.null %>% filter(mean_mille > 0, na.rm=TRUE)
millet.null$logMillet<-log10(millet.null$mean_mille)

#Plot
ggplot(millet.null, aes(x=rescale_ND, y=logMillet)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Millet Near Null Model, T01

Millet Near Null Model, T01 

###  Model - Millet T01
# regular OLS no variance structure
millet.ols <- gls(logMillet ~ rescale_ND, data = millet.null)
# varFixed (variance changes linearly with X)
millet.fixed <- update(millet.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
millet.power <- update(millet.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
millet.exp <- update(millet.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
millet.ConstPower <- update(millet.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(millet.ols, millet.fixed, millet.power, millet.ConstPower, millet.exp)
##                   df      AIC
## millet.ols         3 12808.62
## millet.fixed       3 13255.59
## millet.power       4 12805.30
## millet.ConstPower  5 12807.30
## millet.exp         4 12807.72
#varPower
summary(millet.power)
## Generalized least squares fit by REML
##   Model: logMillet ~ rescale_ND 
##   Data: millet.null 
##       AIC      BIC    logLik
##   12805.3 12832.88 -6398.648
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##     power 
## 0.0553718 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept) -0.16989423 0.026443817 -6.424724       0
## rescale_ND  -0.01371859 0.001853118 -7.402976       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.966
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.8936737 -0.4107742  0.2698067  0.6705882  1.7464964 
## 
## Residual standard error: 0.5033066 
## Degrees of freedom: 7299 total; 7297 residual

#Rapeseed -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
rapeseed.null<-null %>% select(mean_rapes,rescale_ND)
rapeseed.null<-rapeseed.null %>% filter(mean_rapes > 0, na.rm=TRUE)
rapeseed.null$logRapeseed<-log10(rapeseed.null$mean_rapes)

#Plot
ggplot(rapeseed.null, aes(x=rescale_ND, y=logRapeseed)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rapeseed Near Null Model, T01

Rapeseed Near Null Model, T01 

###  Model - Rapeseed T01
# regular OLS no variance structure
rapeseed.ols <- gls(logRapeseed ~ rescale_ND, data = rapeseed.null)
# varFixed (variance changes linearly with X)
rapeseed.fixed <- update(rapeseed.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rapeseed.power <- update(rapeseed.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rapeseed.exp <- update(rapeseed.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rapeseed.ConstPower <- update(rapeseed.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rapeseed.ols, rapeseed.fixed, rapeseed.power, rapeseed.ConstPower, rapeseed.exp)
##                     df      AIC
## rapeseed.ols         3 14872.85
## rapeseed.fixed       3 15691.04
## rapeseed.power       4 14852.02
## rapeseed.ConstPower  5 14853.45
## rapeseed.exp         4 14853.28
#varPower
summary(rapeseed.power)
## Generalized least squares fit by REML
##   Model: logRapeseed ~ rescale_ND 
##   Data: rapeseed.null 
##        AIC      BIC   logLik
##   14852.02 14880.55 -7422.01
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      power 
## 0.07621067 
## 
## Coefficients:
##                   Value   Std.Error    t-value p-value
## (Intercept) -0.09282676 0.015800536  -5.874912       0
## rescale_ND  -0.01535533 0.001219035 -12.596297       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.935
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.2768892 -0.4279157  0.3348087  0.6379341  1.5684906 
## 
## Residual standard error: 0.447725 
## Degrees of freedom: 9259 total; 9257 residual

#Rice -T01

#Select Data and Rescale Distance
null$ric<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
rice.null<-null %>% select(mean_rice_,rescale_ND)
rice.null<-rice.null %>% filter(mean_rice_ > 0, na.rm=TRUE)
rice.null$logRice<-log10(rice.null$mean_rice_)

#Plot
ggplot(rice.null, aes(x=rescale_ND, y=logRice)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rice Near Null Model, T01

Rice Near Null Model, T01 

###  Model - Rice T01
# regular OLS no variance structure
rice.ols <- gls(logRice ~ rescale_ND, data = rice.null)
# varFixed (variance changes linearly with X)
rice.fixed <- update(rice.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rice.power <- update(rice.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rice.exp <- update(rice.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rice.ConstPower <- update(rice.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rice.ols, rice.fixed, rice.power, rice.ConstPower, rice.exp)
##                 df      AIC
## rice.ols         3 13881.81
## rice.fixed       3 15470.57
## rice.power       4 13881.02
## rice.ConstPower  5 13883.02
## rice.exp         4 13880.58
#varPower
summary(rice.power)
## Generalized least squares fit by REML
##   Model: logRice ~ rescale_ND 
##   Data: rice.null 
##        AIC      BIC    logLik
##   13881.02 13909.18 -6936.512
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       power 
## -0.02412247 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept) -0.02183335 0.017118768 -1.275404  0.2022
## rescale_ND   0.01065429 0.001261042  8.448799  0.0000
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.937
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.8420834 -0.4142408  0.2172051  0.6685366  1.6159786 
## 
## Residual standard error: 0.5835948 
## Degrees of freedom: 8428 total; 8426 residual

#Rye -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
rye.null<-null %>% select(mean_rye_Y,rescale_ND)
rye.null<-rye.null %>% filter(mean_rye_Y > 0, na.rm=TRUE)
rye.null$logRye<-log10(rye.null$mean_rye_Y)

#Plot
ggplot(rye.null, aes(x=rescale_ND, y=logRye)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Rye Near Null Model, T01

Rye Near Null Model, T01 

###  Model - Rye T01
# regular OLS no variance structure
rye.ols <- gls(logRye ~ rescale_ND, data = rye.null)
# varFixed (variance changes linearly with X)
rye.fixed <- update(rye.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rye.power <- update(rye.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rye.exp <- update(rye.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rye.ConstPower <- update(rye.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rye.ols, rye.fixed, rye.power, rye.ConstPower, rye.exp)
##                df      AIC
## rye.ols         3 14778.26
## rye.fixed       3 15379.32
## rye.power       4 14756.76
## rye.ConstPower  5 14758.76
## rye.exp         4 14772.60
#varPower
summary(rye.power)
## Generalized least squares fit by REML
##   Model: logRye ~ rescale_ND 
##   Data: rye.null 
##        AIC      BIC    logLik
##   14756.76 14784.91 -7374.381
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      power 
## 0.08782291 
## 
## Coefficients:
##                  Value   Std.Error   t-value p-value
## (Intercept) -0.3955762 0.018380040 -21.52205       0
## rescale_ND   0.0175924 0.001446465  12.16231       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.939
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.7054619 -0.5228168  0.3049936  0.6535288  1.8361489 
## 
## Residual standard error: 0.4701497 
## Degrees of freedom: 8400 total; 8398 residual

#Sorghum -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
sorghum.null<-null %>% select(mean_sorgh,rescale_ND)
sorghum.null<-sorghum.null %>% filter(mean_sorgh > 0, na.rm=TRUE)
sorghum.null$logSorghum<-log10(sorghum.null$mean_sorgh)

#Plot
ggplot(sorghum.null, aes(x=rescale_ND, y=logSorghum)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sorghum Near Null Model, T01

Sorghum Near Null Model, T01 

###  Model - Sorghum T01
# regular OLS no variance structure
sorghum.ols <- gls(logSorghum ~ rescale_ND, data = sorghum.null)
# varFixed (variance changes linearly with X)
sorghum.fixed <- update(sorghum.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sorghum.power <- update(sorghum.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sorghum.exp <- update(sorghum.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sorghum.ConstPower <- update(sorghum.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sorghum.ols, sorghum.fixed, sorghum.power, sorghum.ConstPower, sorghum.exp)
##                    df      AIC
## sorghum.ols         3 18645.01
## sorghum.fixed       3 18952.22
## sorghum.power       4 18477.67
## sorghum.ConstPower  5 18479.67
## sorghum.exp         4 18506.01
#varPower
summary(sorghum.power)
## Generalized least squares fit by REML
##   Model: logSorghum ~ rescale_ND 
##   Data: sorghum.null 
##        AIC     BIC    logLik
##   18477.67 18506.6 -9234.833
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##     power 
## 0.1996702 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept)  0.24963025 0.015136174  16.49229       0
## rescale_ND  -0.03253063 0.001207039 -26.95078       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.923
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.9640782 -0.4787471  0.2296422  0.7081050  1.8790790 
## 
## Residual standard error: 0.3669401 
## Degrees of freedom: 10230 total; 10228 residual

#Soybean -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
soybean.null<-null %>% select(mean_soybe,rescale_ND)
soybean.null<-soybean.null %>% filter(mean_soybe > 0, na.rm=TRUE)
soybean.null$logSoybean<-log10(soybean.null$mean_soybe)

#Plot
ggplot(soybean.null, aes(x=rescale_ND, y=logSoybean)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Soybean Near Null Model, T01

Soybean Near Null Model, T01 

###  Model - Soybean T01
# regular OLS no variance structure
soybean.ols <- gls(logSoybean ~ rescale_ND, data = soybean.null)
# varFixed (variance changes linearly with X)
soybean.fixed <- update(soybean.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
soybean.power <- update(soybean.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
soybean.exp <- update(soybean.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
soybean.ConstPower <- update(soybean.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(soybean.ols, soybean.fixed, soybean.power, soybean.ConstPower, soybean.exp)
##                    df      AIC
## soybean.ols         3 17381.04
## soybean.fixed       3 18239.75
## soybean.power       4 17340.76
## soybean.ConstPower  5 17342.76
## soybean.exp         4 17344.47
#varPower
summary(soybean.power)
## Generalized least squares fit by REML
##   Model: logSoybean ~ rescale_ND 
##   Data: soybean.null 
##        AIC      BIC    logLik
##   17340.76 17369.82 -8666.382
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      power 
## 0.09661652 
## 
## Coefficients:
##                   Value   Std.Error    t-value p-value
## (Intercept)  0.05650602 0.015105583   3.740738   2e-04
## rescale_ND  -0.02214085 0.001149337 -19.264017   0e+00
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.935
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -5.3616282 -0.4015994  0.3470334  0.6821040  1.6531372 
## 
## Residual standard error: 0.4338111 
## Degrees of freedom: 10549 total; 10547 residual

#Sunflower -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
sunflower.null<-null %>% select(mean_sunfl,rescale_ND)
sunflower.null<-sunflower.null %>% filter(mean_sunfl > 0, na.rm=TRUE)
sunflower.null$logSunflower<-log10(sunflower.null$mean_sunfl)

#Plot
ggplot(sunflower.null, aes(x=rescale_ND, y=logSunflower)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Sunflower Near Null Model, T01

Sunflower Near Null Model, T01 

###  Model - Sunflower T01
# regular OLS no variance structure
sunflower.ols <- gls(logSunflower ~ rescale_ND, data = sunflower.null)
# varFixed (variance changes linearly with X)
sunflower.fixed <- update(sunflower.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sunflower.power <- update(sunflower.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sunflower.exp <- update(sunflower.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sunflower.ConstPower <- update(sunflower.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sunflower.ols, sunflower.fixed, sunflower.power, sunflower.ConstPower, sunflower.exp)
##                      df      AIC
## sunflower.ols         3 15781.89
## sunflower.fixed       3 16014.89
## sunflower.power       4 15567.05
## sunflower.ConstPower  5 15545.23
## sunflower.exp         4 15546.02
#varPower
summary(sunflower.power)
## Generalized least squares fit by REML
##   Model: logSunflower ~ rescale_ND 
##   Data: sunflower.null 
##        AIC     BIC    logLik
##   15567.05 15595.6 -7779.525
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##     power 
## 0.2137333 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept) -0.19299200 0.013591393 -14.19957       0
## rescale_ND  -0.01446698 0.001141075 -12.67838       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.907
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.5243282 -0.4099035  0.2536765  0.7261783  1.6669870 
## 
## Residual standard error: 0.3364162 
## Degrees of freedom: 9303 total; 9301 residual

#Wheat -T01

#Select Data and Rescale Distance
null<-read.csv(file="./Near_Null/GCO_Near_T01.csv", header=T)
null$rescale_ND <- null$NEAR_DIST/(1000*1000)
wheat.null<-null %>% select(mean_wheat,rescale_ND)
wheat.null<-wheat.null %>% filter(mean_wheat > 0, na.rm=TRUE)
wheat.null$logWheat<-log10(wheat.null$mean_wheat)

#Plot
ggplot(wheat.null, aes(x=rescale_ND, y=logWheat)) +geom_point() +geom_smooth(method="lm") +theme_justin +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Yield")
## `geom_smooth()` using formula 'y ~ x'
Wheat Near Null Model, T01

Wheat Near Null Model, T01 

###  Model - Wheat T01
# regular OLS no variance structure
wheat.ols <- gls(logWheat ~ rescale_ND, data = wheat.null)
# varFixed (variance changes linearly with X)
wheat.fixed <- update(wheat.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
wheat.power <- update(wheat.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
wheat.exp <- update(wheat.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
wheat.ConstPower <- update(wheat.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(wheat.ols, wheat.fixed, wheat.power, wheat.ConstPower, wheat.exp)
##                  df      AIC
## wheat.ols         3 22286.50
## wheat.fixed       3 23485.16
## wheat.power       4 22266.17
## wheat.ConstPower  5 22267.78
## wheat.exp         4 22267.96
#varPower
summary(wheat.power)
## Generalized least squares fit by REML
##   Model: logWheat ~ rescale_ND 
##   Data: wheat.null 
##        AIC   BIC    logLik
##   22266.17 22296 -11129.09
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      power 
## 0.06416588 
## 
## Coefficients:
##                   Value   Std.Error   t-value p-value
## (Intercept)  0.05869809 0.014548018  4.034783   1e-04
## rescale_ND  -0.00866151 0.001127852 -7.679655   0e+00
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.937
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.7828632 -0.4473925  0.2572152  0.6835164  1.7102002 
## 
## Residual standard error: 0.4933548 
## Degrees of freedom: 12812 total; 12810 residual